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Abstract.  Hysteresis  in  smart  materials  hinders  their  wider  applicabil¬ 
ity  in  actuators.  The  low  dimensional  hysteresis  models  for  these  mate¬ 
rials  are  hybrid  systems  with  both  controlled  switching  and  autonomous 
switching.  In  particular,  they  belong  to  the  class  of  Duhem  hysteresis 
models  and  can  be  formulated  as  systems  with  both  continuous  and 
switching  controls.  In  this  paper,  we  study  the  control  methodology 
for  smart  actuators  through  the  example  of  controlling  a  commercially 
available  magnetostrictive  actuator.  For  illustrative  purposes,  an  infi¬ 
nite  horizon  optimal  control  problem  is  considered.  We  show  that  the 
value  function  satisfies  a  Hamilton-Jacobi-Bellman  equation  (HJB)  of  a 
hybrid  form  in  the  viscosity  sense.  We  further  prove  uniqueness  of  the 
viscosity  solution  to  the  (HJB),  and  provide  a  numerical  scheme  to  ap¬ 
proximate  the  solution  together  with  a  sub-optimal  controller  synthesis 
method.  Numerical  and  experimental  results  based  on  this  approach  are 
presented. 


1  Introduction 

Materials  with  the  intrinsic  characteristics  of  built-in  sensors,  actuators,  and 
control  mechanism  in  their  microstructures,  are  called  smart  materials.  Smart 
materials  and  smart  structures  have  been  receiving  tremendous  interest  in  the 
past  decade,  due  to  their  broad  applications  in  areas  of  aerospace,  manufac¬ 
turing,  defense,  and  civil  infrastructure  systems,  to  name  a  few.  Hysteresis  in 
smart  materials,  e.g.,  magnetostrictives,  piezoceramics,  and  shape  memory  alloys 
(SMAs),  hinders  the  wider  applicability  of  such  materials  in  actuators.  Hystere¬ 
sis  models  can  be  classified  into  phenomenological  models  and  physics-based 
models.  The  most  popular  phenomenological  hysteresis  model  used  in  control 
of  smart  actuators  has  been  the  Preisach  model,  see  e.g.,  [1-3].  An  example  of 
physics-based  model  is  the  Jiles- Atherton  model  for  ferromagnetic  hysteresis  [4], 
where  hysteresis  is  considered  to  arise  from  pinning  of  domain  walls  on  defect 
sites. 

This  paper  is  aimed  at  exploring  the  control  methodology  for  smart  actua¬ 
tors  exhibiting  hysteresis.  We  illustrate  the  ideas  through  the  example  of  con¬ 
trolling  a  commercially  available  magnetostrictive  actuator.  Magnetostriction  is 


the  phenomenon  of  strong  coupling  between  magnetic  properties  and  mechan¬ 
ical  properties  of  some  ferromagnetic  materials  (e.g.,  Terfenol-D):  strains  are 
generated  in  response  to  an  applied  magnetic  field,  while  conversely,  mechani¬ 
cal  stresses  in  the  materials  produce  measurable  changes  in  magnetization.  This 
phenomenon  can  be  used  for  actuation  and  sensing.  Magnetostrictive  actuators 
have  applications  in  micro-positioning,  robotics,  ultrasonics  and  vibration  con¬ 
trol.  Figure  1  shows  a  sectional  view  of  a  Terfenol-D  actuator  manufactured  by 
Etrema  Products,  Inc.  By  varying  the  current  in  the  coil,  we  vary  the  magnetic 
field  in  the  Terfenol-D  rod  and  thus  control  the  motion  of  the  rod  head  (or  the 
force  if  the  motion  is  blocked).  Figure  2  displays  the  hysteresis  observed  in  the 
magnetostrictive  actuator. 


Stainless  Steel  Push  Rod 


Fig.  1.  A  Terfenol-D  actuator  [5]  (Original  source:  Etrema  Products  Inc.) 


Fig.  2.  Hysteresis  in  the  magnetostrictive  actuator 


When  the  frequency  of  the  current  input  I  is  low,  the  magnetostriction  can 
be  related  to  the  bulk  magnetization  M  along  the  rod  direction  through  a  square 
law,  thus  control  of  the  rod  head  boils  down  to  control  of  M.  We  will  employ 
a  low  dimensional  model  [6]  to  represent  the  ferromagnetic  hysteresis  between 


M  and  the  applied  magnetic  field  H ,  where  H  is  proportional  to  /  in  the  low 
frequency  case.  The  model  is  a  hybrid  system  with  both  controlled  switching  and 
autonomous  switching.  It  belongs  to  the  class  of  Duhem  hysteresis  models  and 
can  be  rewritten  as  a  system  involving  both  continuous  control  and  switching 
control.  Note  that  a  low  dimensional  model  for  ferroelectric  hysteresis  has  been 
proposed  by  Smith  and  Horn  [7].  Their  model  has  the  same  structure  as  that  of 
the  ferromagnetic  hysteresis  model  we  use  in  this  paper.  Therefore  the  approach 
we  present  in  this  paper  will  be  fully  applicable  to  control  of  smart  actuators 
made  of  ferroelectric  materials,  e.g.,  piezoelectrics  and  electrostrictives. 

Witsenhausen  formulated  a  class  of  hybrid-state  continuous-time  dynamical 
systems  and  studied  an  optimal  control  problem  back  in  1966  [8].  The  Pontryagin 
Maximum  Principle  or  its  variant  was  used  in  optimal  control  for  hybrid  systems 
in  [9, 10].  By  solving  the  Bellman  inequality,  a  lower  bound  on  the  value  function 
and  an  approximation  to  the  optimal  control  law  were  obtained  in  [11, 12].  Yong 
studied  the  optimal  control  problem  for  a  system  with  continuous,  switching 
and  impulse  controls  in  [13].  With  a  unified  model  for  hybrid  control,  Branicky, 
Borkar  and  Mitter  proposed  a  generalized  quasi- variational  inequalities  (GQVI) 
satisfied  by  the  value  function  [14].  In  this  paper,  we  study  optimal  control  of 
smart  actuators  using  a  viscosity  solutions  approach. 

Dynamic  programming  is  one  of  the  most  important  approaches  in  optimal 
control.  When  the  value  function  of  the  control  problem  is  smooth,  we  can  derive 
the  Hamilton-Jacobi-Bellman  equation  (HJB)  from  the  Dynamic  Programming 
Principle  (DPP),  and  in  many  cases,  solving  the  (HJB)  amounts  to  solving  the 
optimal  control  problem.  The  value  function  however,  in  general,  is  not  smooth 
even  for  smooth  systems,  not  to  mention  for  a  hybrid  system,  like  that  in  our 
model.  Crandall  and  Lions  [15]  introduced  the  notion  of  viscosity  solutions  to 
Hamilton- Jacobi  equations.  This  turned  out  to  be  a  very  useful  concept  for 
optimal  control  since  value  functions  of  many  optimal  control  problems  do  satisfy 
the  (HJB)  in  the  viscosity  sense;  and  under  mild  assumptions,  uniqueness  results 
for  viscosity  solutions  hold  [16].  We  will  explore  this  approach  for  control  of  smart 
actuators.  This  paper  provides  some  flavors  of  this  methodology  by  considering 
an  infinite  horizon  control  problem. 

The  paper  is  organized  as  follows.  In  Sect.  2,  we  present  the  hysteresis  model 
and  examine  its  properties.  In  Sect.  3,  we  formulate  an  optimal  control  prob¬ 
lem  based  on  the  hysteresis  model,  and  show  that  the  value  function  satisfies  a 
hybrid  Hamilton-Jacobi-Bellman  equation  (HJB)  in  the  viscosity  sense.  In  Sect. 
4,  We  prove  that  the  (HJB)  admits  a  unique  solution  in  the  class  of  functions 
to  which  the  value  function  belongs.  We  describe  the  numerical  scheme  to  solve 
the  (HJB)  in  Sect.  5.  This  establishes  the  existence  of  a  solution  to  the  (HJB)  as 
well  as  provides  a  method  for  synthesizing  a  sub-optimal  controller.  Numerical 
and  experimental  results  are  also  reported  in  Sect.  5.  Finally,  conclusions  and 
discussions  are  provided  in  Sect.  6. 

Only  key  proofs  are  provided  in  this  paper  due  to  the  space  limitation.  Other 
proofs  can  be  found  in  [17]. 


2  The  Bulk  Ferromagnetic  Hysteresis  Model 


Jiles  and  Atherton  proposed  a  low  dimensional  model  for  ferromagnetic  hystere¬ 
sis,  based  upon  the  quantification  of  energy  losses  due  to  domain  wall  intersec¬ 
tions  with  inclusions  or  pinning  sites  within  the  material  [4].  A  modification  to 
the  Jiles- Atherton  model  was  made  by  Venkataraman  and  Krishnaprasad  with 
rigorous  use  of  the  energy  balancing  principle  [6].  The  resulting  model,  named 
the  bulk  ferromagnetic  hysteresis  model,  has  a  slightly  different  form  from  the 
Jiles- Atherton  model.  Also  based  on  the  energy  balancing  principle,  they  de¬ 
rived  a  bulk  magnetostrictive  hysteresis  model  [18],  where  high  frequency  effects 
are  considered.  At  low  frequencies,  the  magnetostriction  can  be  related  to  the 
bulk  magnetization  through  a  square  law  [5],  thus  control  of  the  bulk  magneti¬ 
zation  amounts  to  control  of  the  magnetostriction.  In  this  paper,  we  will  restrict 
ourselves  to  the  low  frequency  case  to  highlight  the  methodology  of  hysteresis 
control.  Extension  to  the  high  frequency  case  is  straightforward.  We  now  briefly 
outline  the  bulk  ferromagnetic  hysteresis  model. 

For  an  external  magnetic  field  H  and  a  bulk  magnetization  M,  we  define 
He  =  H  +  olM  to  be  the  effective  field,  where  a  is  a  mean  field  parameter 
representing  inter-domain  coupling.  Through  thermodynamic  considerations,  the 
anhysteretic  magnetization  Man  can  be  expressed  as 

M&n(He)  =  Ms(coth(— )  --£-)  =  Ms£(z)  ,  (1) 

a  He 

where  £(•)  is  the  Langevin  function,  £(z)  =  coth(2)  —  i,  with  2  =  Ms  is  the 
saturation  magnetization  of  the  material  and  a  is  a  parameter  characterizing  the 
shape  of  Man  curve. 

Define 


M  d^(z) 
a  -  acMs— 


h(H,M)  = 


f3(H,M)  = 


cfcMs^§^  -  /x0a(Afan(J?e)  -  M) 
k(a  -  acMs^5l)  +  /z0aa(Man(He)  -  M) 
ckMs ^  +  M0a(Man(ge)  -  M) 
k(a  -  acMs^5l)  -  ^0aa(Man(Ho)  -  M) 


where  c  is  the  reversibility  constant,  /io  is  the  permeability  of  vacuum,  k  is  a 
measure  for  the  average  energy  required  to  break  a  pinning  site.  Note  each  /,  is 
smooth  in  H  and  M. 

The  bulk  ferromagnetic  hysteresis  model  is  as  follows  [6] : 


dM 

~dH 


where  i 


1,  dH  <  0,  M  <  Man(He)  or 
dH  >  0,  M  >  MaJHe) 

2,  dH  <  0,  M  >  MaJHe) 

3,  dH  >  0,  M  <  Man (He) 


If  we  define  a  control  u  =  77 ,  the  model  is  rewritten  as 


(  1,  u  <  0,  M  <  Man(77e )  or 

(H  \  _  (  1  \  ,  I  u  >  0,  M  >  MaJHe) 

\MJ  \  fi(H,  M) )  where  *  \  2,  u  <  0,  M  >  Man(He) 

{  3,  u  >  0,  M  <  Man(He) 


(3) 


Remarks: 


—  Note  that  the  control  u  defined  above  is  different  from  the  physical  current 
7  we  apply  to  the  actuator.  The  current  I  is  related  to  the  state  component 
H  by  a  constant  Cq  (the  coil  factor ):  H  =  cqI.  Therefore  from  the  control  u, 
the  current  we  will  apply  is  I(t)  =  7(0)  +  ^  /0*  u(s)ds. 

—  The  switching  depends  on  both  (the  sign  of)  the  continuous  control  u  and  the 
state  (77,  M),  therefore  the  model  (3)  is  a  hybrid  system  with  both  controlled 
switching  and  autonomous  switching  [14,19]. 

We  can  represent  model  (3)  in  a  more  compact  way.  Letting 


12i={(77,  M)  :  M  <  Man(77e)},  172={(77,  M)  :  M  >  Man(77e)}, 


and  x  =  (77,  M),  we  define 


/+(*) 


if  x  €  1?2 
if  x  £ 


and  /_  (a;) 


if  x  £  fl\ 
if  x  £  b?2 


Since  /,,  1  <  i  <  3,  coincide  on  {(77,  M)  :  M  =  Man(77e)},  f+  and  /_  are 
continuous.  We  define  two  continuous  control  sets 


77+  =  {'«  :  uc  >  u  >  0},  U-  =  {u  :  —uc  <  u  <  0}, 


where  uc  >  0  represents  the  operating  bandwidth  constraint  of  the  actuator 
(recall  u  =  CqI).  To  ease  the  presentation,  we  make  the  dependence  of  switching 
on  u  explicit  by  introducing  a  discrete  control  set  D  =  {1,2}. 

Now  the  model  (3)  can  be  represented  as  a  system  with  both  a  continuous 
control  u  and  a  discrete  mode  (switching)  control  d: 


x  =  f(x,  u,  d)  = 


f+(x)u,  u£U+,  if  <7=1 
f-(x)u ,  uGU-,  if  <7=2 


(4) 


The  (state-dependent)  autonomous  switching  has  now  been  incorporated  into 
the  definitions  of  /+,  /_,  thanks  to  the  nice  structure  of  the  physical  model.  Note 
the  model  (4)  belongs  to  the  category  of  Duhem  hysteresis  model  [20]. 

We  can  prove  that  the  model  enjoys  the  following  properties: 


Proposition  1  (Boundedness  of  fi).  If  the  parameters  satisfy: 


A  acMs 
Ti=a - —  >  0  , 

T2  =  k(a - — )  —  2 n0aaMs  >  0  , 


then  0  <  fi  <  Cf ,  i  —  1,2, 3,  for  some  constant  Cf  >  0. 
Proof.  See  [17]. 


(5) 

(6) 

□ 


Remark  :  Conditions  (5)  and  (6)  are  satisfied  for  typical  parameters.  For  exam¬ 
ple,  taking  the  parameters  identified  in  [5],  a  =  1.9  x  10-4,  a  =  190,  k  =  48 
Gauss,  c  =  0.3,  Ms  =  9.89  x  103  Gauss  and  p,0  =  1,  we  calculate  Ti  =  189.8, 
T2  =  8.40  x  103. 

Proposition  2  (Lipshitz  Continuity).  Functions  f+(x)  and  f-(x)  are  Lip- 
shitz  continuous  with  some  Lipshitz  constant  L,  and  f(x,u,d)  is  Lipshitz  con¬ 
tinuous  with  respect  to  x  with  Lipshitz  constant  Lq  =  Luc. 

Proof.  See  [17].  □ 


3  Optimal  Control:  The  (HJB)  and  Viscosity  Solutions 

We  first  formulate  an  infinite  horizon  optimal  control  problem  for  the  system 
(4).  Define  the  cost  functional  with  an  initial  condition  x  and  a  control  pair 
a(-)  =  {d(-),u(-)j  as 

/»oo 

J(x,  a(-))  =  /  l(x(t),u(t))e~xtdt  ,  (7) 

J  o 

where  the  discount  factor  A  >  0.  Note  the  running  cost  is  defined  to  be 
independent  of  the  switching  control  d,  since  this  makes  sense  in  the  context  of 
smart  actuator  control.  We  require  u(-)  to  be  measurable.  This  together  with 
Proposition  2  guarantees  that  (4)  has  a  unique  solution  x(-)  (the  dependence  of 
x(-)  on  x  and  a(-)  is  suppressed  when  no  confusion  arises). 

The  optimal  control  problem  is  to  find  the  value  function 

V(x)  =  inf  J(x,a(-))  , 

«(■) 

and  if  V(x)  is  achievable,  find  the  optimal  control  «*(•)• 

We  make  the  following  assumptions  about  l(-,  •): 

—  (Ai):  l(x,u)  continuous  in  x  and  u,  l(x,u)  >  0,  \/x,w, 

-  (A2):  | l(xi,u)  -  l(x2,u)\  <  Ci(l  +  |xi|  +  |x2|)|xi  -  x2|,  Vu,  for  some  C\  >  0. 

Note  (A2)  includes  the  case  of  quadratic  cost. 

We  can  show  the  value  function  is  locally  bounded  and  locally  Lipshitz  con¬ 
tinuous. 


Proposition  3  (Local  Boundedness).  Under  assumptions  (Hi)  and  (H2), 
VA  >  0,  V(x)  is  locally  bounded,  i.e.,  Vi?  >0,3  Cr  >  0,  such  that  |T> (cc) |  < 

Cr,  Vx  £  B(0,  R)  =  {x  :  |x|  <  R} . 

Proof.  See  [17].  □ 

Proposition  4  (Local  Lipshitz  Continuity).  Under  assumptions  (H-| )  and 
(H2),  VA  >  2Lq  with  Lq  as  defined  in  Proposition  2,  V{x)  is  locally  Lipshitz,  i.e., 
Vi?  >  0,3Lr  >  0,  such  that  \V(x\)  —  V{x2)\  <  Lr|xi  —  X2I,  Vxi,#2  £  B(0,R). 
In  addition,  Lr  can  be  chosen  to  be  C(1  +  R)  for  some  C  >  0. 

Proof.  See  [17].  □ 

Remarks:  The  proof  of  Proposition  4  uses  the  bounds  for  \x\{t)  —  X2(t)\  and 
|xi(t)|,  where  Xi(-),X2(-)  are  two  trajectories  starting  from  aq  and  X2 ■  The 
bounds  are  obtained  using  Proposition  2  and  the  Gronwall  inequality.  One  can 
get  a  sharper  estimate  for  |xi(t)|  (linear  growth)  by  exploiting  Proposition  1. 
This  can  be  used  to  weaken  the  condition  A  >  2Lq  to  A  >  Lq  in  Proposition  4 
and  anywhere  else  it  appears. 

The  value  function  satisfies  the  Dynamic  Programming  Principle  (DPP)  : 
Proposition  5  (DPP).  Assume  (Ai)  and  (A2),  A  >  0.  We  have 

V(x)  =  inf{  f  e~Xsl(x(s),  u(s))ds  +  e~xtV(x(t))},  Vt  >  0,  Vx  .  (8) 

Proof.  The  argument  is  standard,  see  [17].  □ 

Based  on  the  (DPP),  we  can  show  that  the  value  function  P(-)  satisfies 
a  Hamilton- Jacobi-Bellman  equation  (HJB)  of  a  hybrid  type  in  the  viscosity 
sense.  Viscosity  solutions  to  Hamilton- Jacobi  equations  were  first  introduced  by 
Crandall  and  Lions  [15].  Here  we  use  one  of  the  three  equivalent  definitions  [21]: 

Definition  1  (Viscosity  Solution).  Let  W  be  a  continuous  function  from  an 
open  set  O  £  M"  into  R  and  let  DW  denote  the  gradient  of  W  (when  W  is 
differentiable) .  We  call  W  a  viscosity  solution  to  a  nonlinear  first  order  par¬ 
tial  differential  equation  F(x,W(x),DW(x))=  0,  provided  it  is  both  a  viscosity 
subsolution  and  viscosity  super  solution;  and  by  viscosity  sub  (super)  solution,  we 
mean:  V</>  £  C1(0),  if  W  —  tj>  attains  a  local  maximum  (minimum)  at  Xq  £  O, 
then  F[xq,W{xq),  D(f[xf)))<  (>)0. 

Theorem  1  (HJB).  Assume  (Hi)  and  (A2),  A  >  2Lq.  V(x)  is  a  viscosity 
solution  of: 

A W (x)  +  max{  max  {— u/+(x)  •  DW (x)  —  l(x,  it)}, 
ueu+ 

ma x{—uf-(x)  ■  DW(x)  —  l(x,u)}}  =  0  .  (9) 

U£U- 


Proof.  See  [17]. 


□ 


4  Uniqueness  of  the  Solution  to  the  (HJB) 


We  would  like  to  characterize  the  value  function  V  as  a  unique  solution  to  the 
(HJB).  The  uniqueness  result  basically  follows  from  Theorem  1.5  in  [22].  In  [22], 
the  author  gave  only  a  sketch  of  proof.  For  completeness,  we  will  provide  the  full 
proof  here. 

Before  stating  the  theorem,  we  first  identify  structural  properties  of  the 
(HJB).  We  rewrite  (9)  as: 

A W(x)  +  H(x,  DW(x))  =  0  ,  (10) 


where 

H(x,p)  =  max{max{-«/+(a;)  •  p  —  l(x,u)},  max{— uf-(x)  ■  p  —  l(x,u)}}  . 
u€lU- I-  u£U- 


Proposition  6.  Assume  (A2).  H(x,p)  satisfies  the  following: 

\H(x!,p)  -  H(x2,p)\  <Cr(1+\p\)\xi  -  x2\,  Vx!,x2  €  B(Q,R),  Vp  ,  (11) 
\H(x,p1)-H(x,p2)\<C0\pi-p2\,Vx,Vp1,p2  ,  (12) 

for  some  Cr  >  0,Co  >  0,  with  Cr  dependent  on  R. 

Proof.  We  will  only  prove  (11),  since  proof  of  (12)  is  analogous. 

Without  loss  of  generality,  suppose  ui  £  U-  attains  the  maximum  in  H (xi ,  p) . 
Since  H(x2,p)  >  —U\f-(x2)  -p  —  l(x2,u i), 

H(xi,p)  -  H(x2,p)  <  —uif-(xi)  -p  -  l(xi,ui)  +  uif-(x2)  -p  +  l(x2,ui) 

<  \p\L0\xi  -  x2\  +  Ci(l  +  |xi|  +  |ar2|)|a:i  -  x2| 

<  Cr(1  +  \p\)\x i  -  x2\  , 

where  Cr  is  a  constant  dependent  on  R.  By  symmetry,  we  conclude.  □ 


Remarks:  As  we  have  seen  above,  despite  the  hybrid  structure  of  our  physical 
model,  H (x,  p)  enjoys  nice  structural  properties,  which  enables  us  to  prove  the 
uniqueness  result. 

From  Proposition  4,  we  know  that  the  value  function  P(-)  belongs  to  the 
class 


P(K2)  =  {W(-)  :  \W{x\)  —  W{x2)\  <  C(1  +  i?)|xi  —  x2\,Wxi,x2  £  B(0,R), 
for  some  C  >  0}. 

The  following  theorem  is  adapted  from  Theorem  1.5  in  [22]. 

Theorem  2.  If  (10)  has  a  viscosity  solution  in  'P(R2),  it  is  unique. 


Proof.  Without  loss  of  generality,  we  take  A  =  1.  Let  eT^R2)  be 

viscosity  solutions  to  (10).  For  e  >  0,  a  >  0,  m  >  2,  define 

<P(x,y)  =  W(x)  -  V(y)  -  ^  ~  -  a(<  x  >m  +  <  y  >m)  , 

with  <  x  >=  sj\  +  \x\ 2.  Since  gT^R2),  lim|x|+|y|_>00  <P(x,  y)=  — oo. 

By  continuity  of  <P(-,  •),  there  exists  (xq,  yo)  where  <P  attains  the  global  maximum. 
First  we  want  to  obtain  bounds  for  |#o|,  |j/o|  and  \xo  —  2/0 1- 

From  ^(0,0)  <  <P(xo,yo),  and  W{-),V{-)  gT-^R2),  we  can  get 

<  x0  >m  +  <  yo  >m<  Ca{  1+  <x0>2  +  <y0  >2)  , 

where  Ca  is  a  constant  independent  of  e  (but  dependent  on  a).  Since  m  >  2, 
there  exists  Ra  >  0  (independent  of  e),  such  that  |xo|  <  Ra,  |j/o|  <  Ra- 
From  4>(x0,  x0)  +  $(yo,  yo)  <  2  #(x0,  yo),  we  can  derive 

\xo  -yo\<eC'a  ,  (13) 


with  C'a  depending  on  a  only. 

Define 

</>(%)  =  V(y0)  +  ~  \x  -  yo  I2  +  a(<  x  >m  +  <  y0  >m )  , 

Ml /)  =  W(x  0)  -  i|x0  -  j/| 2  -  a(<  x0>m  +  <y  >m )  . 

Since  W  —  (f>  achieves  maximum  at  Xq,  and  V  —  ip  achieves  minimum  at  yo, 

W(xo)  +  H(xo,Dcf(xo))<0  ,  (14) 

V(y0)  +  H(y0,Dip(yo))  >  0  .  (15) 

Subtracting  (15)  from  (14)  and  using  Proposition  6,  we  have 

W(x0)  -V(y0)  <  CRa(l  +  ~  |x0  -  yo  1)1*0  -  Vo\ 

+aC0m{<  xo  >m_1  +  <  2/0  >m_1)  • 

Now  fix  a,  construct  a  sequence  {efc }  with  lim^oo  efc  =  0.  We  denote  the 
corresponding  maximizers  of  $  as  ( Xok,Vok )•  Since  \/k,  (xok,  2/ofc)  G  B(0,Ra),  by 
extracting  a  subsequence  if  necessary,  we  get 

lim  (xofe,  2/ofc)  — ►  (*a,  2/a)  e  B(0,Ra)  ■  (16) 

k—>oo 

From  (13),  we  have  xa  =  ya.  For  each  efc,  from  <P(x ,  x)  <  @(xok,yok),  we  can  get 

2 

W(x)  -  V{x)  -2  a<x  >m<  CRa(  1  H - |x0fc  -  yok\)\x0k  -  2/ofc  | 

Cfc 

+  aC0m(<  x0 k  >m~ 1  +  <  yok  >m_1)  -  a(<  x0k  >m  +  <  yok  >m )  , 


and  letting  k  — »  oo, 

W{ x)  —  V(x)  <  2a(C0m  <  xa  >m~1  —  <  xa  >m )  +  2a  <  x  >m  . 

Since  Com  <  xa  >m_1  —  <  xa  >m<  C"  for  some  C"  >  0, 

W{x)  -  V(x)  <  2 a{C"+  <  x  >m)  . 

Letting  a  — >  0,  we  get  W(x)  —  V(x)  <  0,  \/x.  We  conclude  by  noting  W  and  V 
are  symmetric.  □ 

From  Theorem  2,  if  we  can  solve  for  a  solution  to  (10)  in  ^(K2),  it  must  be  the 
value  function.  One  way  to  solve  it  is  by  discrete-time  approximation. 

5  The  Discrete  Approximation  Scheme 

The  approximation  will  be  accomplished  in  two  steps.  First  we  approximate  the 
continuous  time  optimal  control  problem  by  a  discrete  time  problem,  derive  the 
hybrid  discrete  Bellman  equation  (DBE),  and  show  the  value  function  of  the 
discrete  time  problem  converges  to  that  of  the  continuous  time  problem  locally 
uniformly.  Following  [16],  we  call  this  step  “semi-discrete”  approximation.  Then 
we  indicate  how  to  further  discretize  (DBE)  in  the  spatial  variable,  which  is  called 
“fully-discrete”  approximation.  The  approaches  we  take  here  follow  closely  those 
in  [16] (Chapter  VI  and  Appendix  A). 

Consider  a  discrete  time  problem  obtained  by  discretizing  the  original  con¬ 
tinuous  time  one  with  time  step  h  £  (0,  j).  The  dynamics  is  given  by 

x[n]  =  x[n  —  1]  +  hf(x[n  —  1],  u[n  —  1],  d[n  —  1]),  ar[0]  =  x  ,  (17) 

and  the  cost  is  given  by 

OO 

Jh{x,a[])  =  Y  hl(x[n],u[n])(  1  -  A h)n  ,  (18) 

71= 0 

where  «[■]  =  {d[-],it[-]}  is  the  control.  The  value  function  is  defined  to  be 

Vh(x)  =  inf  Jh(x,  «[•])  .  (19) 

«[•] 

It’s  not  hard  to  show: 

Proposition  7.  Assume  A±  and  A2,  A  >  2L0.  Then  I4(-)  €  'P(M2),  and  the 
coefficient  C  in  defining  P(R2)  can  be  made  independent  of  h. 

Following  standard  arguments,  one  can  show: 

Proposition  8  (DBE).  V),, ( ■ )  satisfies: 

Vh(x)  =  min{  min  {(1  —  Xh)Vh(x  +  huf+{x))  +  hl( x,  u)}, 

min  {(1  —  Xh)Vh(x  +  huf-(x))  +  hl(x,  w)}}  .  (20) 

uGU- 


It’s  of  interest  to  know  whether  (20)  characterizes  the  value  function  Vh(x). 
Unlike  in  [16](Chapter  VI),  where  a  bounded  value  function  was  considered,  we 
have  Vh  unbounded.  But  it  turns  out  that  with  a  little  bit  additional  assumption, 
(20)  has  a  unique  solution. 

Proposition  9.  There  exists  a  unique  solution  in  "P(R2)  to  (20),  if 


Q.-Xh)(y/CffZ+C0) 

v'cf—\  Co 


where  Cq  =  huc 


1 

C{ 


and  Cf  is  as  defined  in  Proposition  1. 


(21) 


Proof.  Let  14 (x)  =  Vh{x)  <  x  >  m,ra  >  2,  where  <  x  >=  yU  +  |cc|2.  Since 
Vh  €  P(R2),  Vh  is  bounded.  In  terms  of  14,  (20)  is  rewritten  as 

Vh(x)  =  (g(Vh))(x)=mm{  (22) 

min  {(1  —  \h)Vh(x  +  huf+(x))  <  1  +  ^uf+(x)  > - ^  mix, u\  <  x  >_m}, 

uGU+  <  x  >m 

min  {(1  —  \h))Vh{x  +  kuf-{x))  <  —  > - 1_  <  x  >-m}}  . 

ueU-  <  x  >m 

It  suffices  to  show  (22)  has  a  unique  solution.  It’s  clear  that  the  operator  £/(•) 
maps  any  W  €  BC( R2)  into  BC{ R2),  where  I?C(R2)  denotes  the  set  of  bounded 
continuous  functions.  When  (21)  is  satisfied,  one  can  show  that  G(-)  is  a  contrac¬ 
tion  mapping.  Hence  we  conclude  using  the  Contraction  Mapping  Principle.  □ 

The  following  theorem  asserts  that  I4(-)  converges  to  U(-)  as  h  — >  0.  The 
proof  can  be  found  in  [16] (Chapter  VI) (with  minor  modification). 

Theorem  3.  Assume  Ai  and  A2,  A  >  2L0,  and  (21).  Then 

sup  \Vh(x)  —  V(a:)|  — ►  0  as  h  —>  0,  (23) 

xeK, 

for  every  compact  K, C  R2. 

It  was  also  shown  in  [16]  that  one  can  obtain  a  sub-optimal  control  for  the 
continuous  time  problem  when  solving  the  (DBE) .  Theoretically  the  solution  to 
(20)  can  be  obtained  by  successive  approximation.  A  practical  approximation 
scheme  for  solving  the  (DBE)  is  described  in  [16] (Appendix  A,  by  Falcone), 
which  we  have  followed  in  the  numerical  simulation.  It  was  shown  there  that 
when  space  discretization  gets  finer  and  finer,  the  solution  obtained  via  solving 
a  finite  system  of  equations  converges  to  I4(-). 

Computation  can  only  be  done  in  a  bounded  domain.  The  domain  we  used 
in  simulation  is  of  the  form  17=  {( H,M )  :Hm ;n<  H  <  I7max,  |M|  <  Ms}.  The 
constraint  \M\  <  Ms  arises  from  the  physics,  while  the  constraint  on  H  is  due 


to  limitation  on  the  range  of  current  input.  The  value  function  of  an  optimal 
control  problem  with  state-space  constraints  is  a  constrained  viscosity  solution  to 
the  (HJB)  in  ft  [23] ,  namely,  a  solution  in  the  interior  of  ft  and  a  supersolution 
in  ft.  Theoretical  results  for  the  constrained  state-space  case  are  omitted  in  this 
paper. 

The  values  we  used  for  the  model  parameters  are  as  those  in  the  remarks 
following  Proposition  1.  The  running  cost  was  defined  as  l(H,  M,u)=100(H  — 
Hq)2+0.1M2+0.01u2,  where  Hq  corresponds  to  some  desired  steady  current 
input.  Other  parameters:  Hmm  =  19.8,  Hmax  =  407.8,  H0  =  213.8,  A  =  1.58  x 
103,  h  =  5  x  10-4,  uc  =1.22  x  103.  Each  of  £7+  and  U-  is  discretized  into  20 
levels,  while  each  dimension  of  the  state  space  is  discretized  into  40  levels. 

Figure  3  shows  the  value  function  and  the  optimal  feedback  control  map. 
Optimal  trajectories  obtained  through  simulation  and  experiments  from  three 
different  initial  conditions  (A,  B,  C)  are  shown  in  Fig.  4,  where  the  arrows 
indicate  directions  of  evolution  as  well  as  stationary  points  of  the  closed-loop 
systems.  Figure  5  shows  the  experimental  setup. 


Fig.  3.  (a)  Value  function  and  its  level  curves  (b)  Optimal  feedback  control  map 


6  Conclusions  and  Discussions 

In  this  paper,  we  have  studied  the  viscosity  solutions  approach  for  optimal  control 
of  a  smart  actuator  based  on  the  low  dimensional  hysteresis  model.  We  took  the 
infinite  horizon  optimal  control  problem  as  an  example  and  characterized  the 
value  function  as  the  (unique)  solution  of  a  Hamilton- Jacobi-Bellman  equation 
of  a  hybrid  form.  We  pointed  out  how  to  solve  the  (HJB)  and  obtain  a  sub- 
optimal  control  by  discrete  time  approximation. 

Future  work  includes  extension  of  this  approach  to  other  control  problems  of 
practical  interest,  including  the  finite  horizon  control  problem,  the  time-optimal 
control  problem  and  the  H0 0  control  problem.  Also  the  stability  issues  associated 
with  the  closed-loop  systems  need  to  be  investigated. 


H(Oe) 


Fig.  4.  Optimal  trajectories 


Amplifier 


Fig.  5.  Experimental  setup 
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